import matplotlib.pyplot as plt
from numpy import concatenate
from numpy import flipud
from numpy import sum
from numpy import zeros
from numpy.random import rand
from pywt import Wavelet
from pywt import dwt
from scipy.signal import cont2discrete
from scipy.signal import impulse
from scipy.signal import lfilter

u = rand(100)

u1 = u[1:68]
last = len(u1);

trans = 'db4';

[lo_d, hi_d, lo_r, hi_r] = Wavelet(trans).filter_bank
nc = len(lo_d)-1
den = concatenate(([1], zeros(nc)), axis=1)

[ca1, cd1] = dwt(u1, trans)

U1 = concatenate((u1, zeros(nc)), axis=0)

tot = len(U1)

y1 = lfilter(lo_d, den, U1)
y2 = lfilter(hi_d, den, U1)

cca1 = y1[2:2:tot]
ccd1 = y2[2:2:tot]

print lo_d
print den
print impulse((lo_d, den), N=tot)
h1 = impulse((lo_d, den), N=tot)

fig = plt.figure(num=10, figsize=(9, 14), dpi=80, facecolor='w', edgecolor='k')
fig.suptitle("frequency response of the filters")
plt.subplot(111)
plt.title(r'Decomposition LP filter: H')
plt.plot(h1)
plt.show()

h2 = impulse(cont2discrete(hi_d, den), n=tot)

for i in range(30):
    u2 = u[i + 1:i + 68]
    U2 = [u2, zeros(nc, 1)]

    y3 = zeros(tot, 1)
    y4 = y3

    y3[1:last-1] = y1[2:last]-h1[2:last] * u1[1]
    y3[last:tot-1] = y1[last + 1:tot]-h1[last + 1] * u1[1] + h1[1:tot-last] * u2[last]
    y3[tot] = sum(h1 * flipud(U2))


    y4[1:last-1] = y2[2:last]-h2[2:last] * u1[1]
    y4[last] = sum(h2 * flipud(U2));

    [ca2, cd2] = dwt(u2, trans)

    cca2 = y3[2:2:tot]
    ccd2 = y4[2:2:tot]

    cca2 = cca2[:len(ca2)]

    u1 = u2
    y1 = y3
    y2 = y4
